Гипотеза: если полигоны закрытой тесселяции зданий (т.е. ограниченные дорожной сетью) соприкасаются с дорожной сетью уровня motorway, trunk, primary, secondary, tertiary, то значит здание находится на первой линии, если нет соприкосновения - здание находится на второй линии.
# ! pip install shapely
# ! pip install osmnx
# ! pip install momepy
# ! pip install geopandas
# ! pip install pandas
# ! pip install numpy
# ! pip install matplotlib
# ! pip install folium
# ! pip install mapclassify
# ! pip install shapely==1.8.5.post1
# ! pip install osmnx==1.3.0
# ! pip install momepy==0.6.0
# ! pip install geopandas==0.12.2
# ! pip install pandas==1.3.4
# ! pip install numpy==1.20.3
# ! pip install matplotlib
# ! pip install folium
# ! pip install mapclassify
# импорт библиотек
import osmnx as ox
import momepy as mm
import geopandas as gpd
import pandas as pd
import numpy as np
import folium
from mapclassify import classify
from matplotlib import pyplot as plt
import shapely
from shapely.geometry import Polygon, LineString, MultiPoint, Point
ox.__version__
'1.9.1'
shapely.__version__
'2.0.2'
mm.__version__
'0.7.0'
# загрузка из OSM зданий внутри границ города
gdf = ox.geometries.geometries_from_place('Kahla, Germany', tags={'building':True}) #
C:\Users\Asus\AppData\Local\Temp\ipykernel_18436\2177612656.py:2: FutureWarning: The `geometries` module and `geometries_from_X` functions have been renamed the `features` module and `features_from_X` functions. Use these instead. The `geometries` module and function names are deprecated and will be removed in the v2.0.0 release.
gdf = ox.geometries.geometries_from_place('Kahla, Germany', tags={'building':True}) #
# проецирование в UTM зону геодатафрейма со зданиями
gdf_projected = ox.projection.project_gdf(gdf)
gdf_projected = gdf_projected[~gdf_projected.geometry.is_empty]
gdf_projected.plot()
<Axes: >
gdf_projected.index
MultiIndex([('way', 45059580),
('way', 87337221),
('way', 87337222),
('way', 87337224),
('way', 87337225),
('way', 87337226),
('way', 87337230),
('way', 87337233),
('way', 87339964),
('way', 87371536),
...
('way', 1110828827),
('way', 1110828828),
('way', 1110828829),
('way', 1110828830),
('way', 1110828831),
('way', 1110828832),
('way', 1112005241),
('way', 1112005242),
('way', 1112005243),
('way', 1112005244)],
names=['element_type', 'osmid'], length=3013)
# подготовка полигонов со зданиями для дальнейшей работы
buildings = mm.preprocess(gdf_projected, size=30,compactness=0.2, islands=True, verbose=False)
C:\Users\Asus\anaconda3\envs\geo\Lib\site-packages\momepy\preprocessing.py:84: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning sw = libpysal.weights.contiguity.Rook.from_dataframe(blg, silence_warnings=True) C:\Users\Asus\anaconda3\envs\geo\Lib\site-packages\momepy\preprocessing.py:84: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning sw = libpysal.weights.contiguity.Rook.from_dataframe(blg, silence_warnings=True)
# загрузка из OSM дорожной сети внутри границ города
streets_graph = ox.graph_from_place('Kahla, Germany', network_type='drive') # Kahla, Germany
# проецирование в UTM зону геодатафрейма с дорожной сетью
streets_graph = ox.projection.project_graph(streets_graph)
# создание геодатафрейма дорожной сети из графа
streets = ox.graph_to_gdfs(streets_graph, nodes=False, edges=True, node_geometry=False, fill_edge_geometry=True)
# визуализация дорожной сети и зданий
ax = buildings.plot(figsize=(10, 10))
streets.plot(ax=ax, color='k', linewidth=.5)
ax.set_axis_off()
# задаем ограничение для прорисовки полигонов тесселяции зданий
limit = mm.buffered_limit(buildings)
# создаем полигоны пересечений улиц
enclosures = mm.enclosures(streets, limit=gpd.GeoSeries([limit]))
C:\Users\Asus\anaconda3\envs\geo\Lib\site-packages\geopandas\array.py:1459: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as WGS 84 / UTM zone 32N (the single non-null crs provided). return GeometryArray(data, crs=_get_common_crs(to_concat))
enclosures.plot()
<Axes: >
df_streets = streets.reset_index()
df_streets
| u | v | key | osmid | name | highway | maxspeed | oneway | reversed | length | lanes | ref | geometry | bridge | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 31918054 | 313327901 | 0 | 154934607 | Gerberstraße | residential | 20 | False | False | 6.545 | NaN | NaN | LINESTRING (681804.922 5630396.235, 681806.646... | NaN |
| 1 | 31918054 | 1583262490 | 0 | [370409322, 370409324] | Rudolstädter Straße | primary | 50 | False | True | 456.417 | [3, 2] | B 88 | LINESTRING (681804.922 5630396.235, 681840.863... | NaN |
| 2 | 313327901 | 31918054 | 0 | 154934607 | Gerberstraße | residential | 20 | False | True | 6.545 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681804.922... | NaN |
| 3 | 313327901 | 1673865104 | 0 | 154934607 | Gerberstraße | residential | 20 | False | False | 679.145 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681815.533... | NaN |
| 4 | 313327901 | 313328313 | 0 | [154934608, 363472734] | Rudolstädter Straße | residential | NaN | False | False | 70.289 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681806.472... | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 390 | 2179857603 | 2179857601 | 0 | 207735308 | Am Aschborn | residential | NaN | False | False | 32.753 | NaN | NaN | LINESTRING (681373.057 5632112.685, 681341.840... | NaN |
| 391 | 2179857599 | 2179857594 | 0 | 207735312 | Am Aschborn | residential | NaN | False | True | 29.405 | NaN | NaN | LINESTRING (681414.507 5632080.613, 681386.452... | NaN |
| 392 | 2179857601 | 2179857603 | 0 | 207735308 | Am Aschborn | residential | NaN | False | True | 32.753 | NaN | NaN | LINESTRING (681341.840 5632102.458, 681373.057... | NaN |
| 393 | 2179857606 | 2179857603 | 0 | 154934582 | Am Aschborn | residential | NaN | False | False | 28.886 | NaN | NaN | LINESTRING (681400.427 5632122.176, 681373.057... | NaN |
| 394 | 2195363605 | 2195363607 | 0 | 209309842 | Am Storchenheim | residential | 50 | False | True | 23.521 | NaN | NaN | LINESTRING (682643.098 5630212.669, 682644.372... | NaN |
395 rows × 14 columns
buildings['uID'] = mm.unique_id(buildings)
enclosed_tess = mm.Tessellation(buildings, unique_id='uID',enclosures=enclosures).tessellation
C:\Users\Asus\anaconda3\envs\geo\Lib\site-packages\momepy\elements.py:451: FutureWarning: The `query_bulk()` method is deprecated and will be removed in GeoPandas 1.0. You can use the `query()` method instead. inp, res = buildings.sindex.query_bulk( C:\Users\Asus\AppData\Local\Temp\ipykernel_18436\3524002744.py:1: UserWarning: dask.dataframe could not be imported. Setting `use_dask=False`. enclosed_tess = mm.Tessellation(buildings, unique_id='uID',enclosures=enclosures).tessellation
enclosed_tess.shape
(2580, 3)
ax = enclosed_tess.plot(figsize=(10, 10))
buildings.plot(ax=ax, color='white', alpha=.5)
ax.set_axis_off()
# визуализируем границы enclosed tesselations и enclosures (красная линия)
ax = enclosed_tess.plot(edgecolor='white', linewidth=0.2, figsize=(10, 10))
buildings.plot(ax=ax, color='white', alpha=.5)
enclosures.plot(ax=ax, edgecolor='red', linewidth=2, alpha=.2)
ax.set_axis_off()
ax.set_xlim(681500, 682500)
ax.set_ylim(5631000, 5632000)
(5631000.0, 5632000.0)
enclosed_tess
| uID | geometry | eID | |
|---|---|---|---|
| 0 | 1373.0 | POLYGON ((681928.695 5630536.070, 681928.721 5... | 0 |
| 1 | 1370.0 | POLYGON ((681949.898 5630592.392, 681950.869 5... | 0 |
| 2 | 1365.0 | POLYGON ((682027.202 5630639.526, 682027.100 5... | 0 |
| 3 | 1371.0 | POLYGON ((681930.699 5630535.484, 681929.697 5... | 0 |
| 4 | 1368.0 | POLYGON ((681964.672 5630620.573, 681965.876 5... | 0 |
| ... | ... | ... | ... |
| 2575 | NaN | POLYGON ((681852.809 5632875.693, 681853.585 5... | 56 |
| 2576 | NaN | POLYGON ((682369.443 5630544.364, 682378.978 5... | 57 |
| 2577 | 2153.0 | POLYGON ((680458.914 5630872.499, 680464.377 5... | 58 |
| 2578 | 568.0 | POLYGON ((681002.891 5632199.397, 681000.911 5... | 59 |
| 2579 | 320.0 | POLYGON ((683485.762 5631278.394, 683479.331 5... | 61 |
2580 rows × 3 columns
# для каждого полигона enclosed_tess найдем соприкосновение с улицами и запишем в датафрейм
lst = []
for index, row in enclosed_tess.iterrows():
streets_selected = row['geometry'].touches(df_streets['geometry'])
streets_selected_overlaps = row['geometry'].overlaps(df_streets['geometry'].buffer(1))
df_streets_selected = list(set(list(df_streets[streets_selected].index) + list(df_streets[streets_selected_overlaps].index)))
enclosure_tess_index = index
street_index = df_streets_selected
l = [enclosure_tess_index, street_index]
lst.append(l)
df_enclosure_tess_index = pd.DataFrame(lst, columns = ['enclosure_tess_index', 'street_indexes'])
df_enclosure_tess_index.street_indexes
0 [0, 129, 2, 3, 4, 1, 7, 50]
1 [1, 50, 3, 7]
2 [50, 3]
3 [1, 50, 3, 7]
4 [1, 50, 3, 7]
...
2575 [290, 306]
2576 []
2577 []
2578 []
2579 []
Name: street_indexes, Length: 2580, dtype: object
df_streets
| u | v | key | osmid | name | highway | maxspeed | oneway | reversed | length | lanes | ref | geometry | bridge | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 31918054 | 313327901 | 0 | 154934607 | Gerberstraße | residential | 20 | False | False | 6.545 | NaN | NaN | LINESTRING (681804.922 5630396.235, 681806.646... | NaN |
| 1 | 31918054 | 1583262490 | 0 | [370409322, 370409324] | Rudolstädter Straße | primary | 50 | False | True | 456.417 | [3, 2] | B 88 | LINESTRING (681804.922 5630396.235, 681840.863... | NaN |
| 2 | 313327901 | 31918054 | 0 | 154934607 | Gerberstraße | residential | 20 | False | True | 6.545 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681804.922... | NaN |
| 3 | 313327901 | 1673865104 | 0 | 154934607 | Gerberstraße | residential | 20 | False | False | 679.145 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681815.533... | NaN |
| 4 | 313327901 | 313328313 | 0 | [154934608, 363472734] | Rudolstädter Straße | residential | NaN | False | False | 70.289 | NaN | NaN | LINESTRING (681806.646 5630389.917, 681806.472... | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 390 | 2179857603 | 2179857601 | 0 | 207735308 | Am Aschborn | residential | NaN | False | False | 32.753 | NaN | NaN | LINESTRING (681373.057 5632112.685, 681341.840... | NaN |
| 391 | 2179857599 | 2179857594 | 0 | 207735312 | Am Aschborn | residential | NaN | False | True | 29.405 | NaN | NaN | LINESTRING (681414.507 5632080.613, 681386.452... | NaN |
| 392 | 2179857601 | 2179857603 | 0 | 207735308 | Am Aschborn | residential | NaN | False | True | 32.753 | NaN | NaN | LINESTRING (681341.840 5632102.458, 681373.057... | NaN |
| 393 | 2179857606 | 2179857603 | 0 | 154934582 | Am Aschborn | residential | NaN | False | False | 28.886 | NaN | NaN | LINESTRING (681400.427 5632122.176, 681373.057... | NaN |
| 394 | 2195363605 | 2195363607 | 0 | 209309842 | Am Storchenheim | residential | 50 | False | True | 23.521 | NaN | NaN | LINESTRING (682643.098 5630212.669, 682644.372... | NaN |
395 rows × 14 columns
df_enclosure_tess_index_result = df_enclosure_tess_index.explode('street_indexes')
# типы дорожной сети, принадлежность к которой будет считаться 1й линией домов
first_line_street_types = ['motorway', 'trunk', 'primary', 'secondary', 'tertiary']
streets_eti = pd.merge(df_enclosure_tess_index_result, df_streets[['highway']], right_index=True, left_on='street_indexes', how='left')
streets_eti['is_first_line'] = streets_eti['highway'].map(lambda x: 1 if x in first_line_street_types else 0)
# полигоны enclosed_tess с признаком: 1 линия или нет
streets_eti_prcd = streets_eti.groupby('enclosure_tess_index').agg({'is_first_line': 'max'}).reset_index(drop=False)
streets_eti_prcd.describe()
| enclosure_tess_index | is_first_line | |
|---|---|---|
| count | 2580.000000 | 2580.000000 |
| mean | 1289.500000 | 0.151938 |
| std | 744.926171 | 0.359030 |
| min | 0.000000 | 0.000000 |
| 25% | 644.750000 | 0.000000 |
| 50% | 1289.500000 | 0.000000 |
| 75% | 1934.250000 | 0.000000 |
| max | 2579.000000 | 1.000000 |
streets_eti_prcd.groupby('is_first_line').size()
is_first_line 0 2188 1 392 dtype: int64
df_streets.groupby('highway').size()
highway primary 20 primary_link 1 residential 255 secondary 28 secondary_link 1 tertiary 69 unclassified 21 dtype: int64
enclosed_tess
| uID | geometry | eID | |
|---|---|---|---|
| 0 | 1373.0 | POLYGON ((681928.695 5630536.070, 681928.721 5... | 0 |
| 1 | 1370.0 | POLYGON ((681949.898 5630592.392, 681950.869 5... | 0 |
| 2 | 1365.0 | POLYGON ((682027.202 5630639.526, 682027.100 5... | 0 |
| 3 | 1371.0 | POLYGON ((681930.699 5630535.484, 681929.697 5... | 0 |
| 4 | 1368.0 | POLYGON ((681964.672 5630620.573, 681965.876 5... | 0 |
| ... | ... | ... | ... |
| 2575 | NaN | POLYGON ((681852.809 5632875.693, 681853.585 5... | 56 |
| 2576 | NaN | POLYGON ((682369.443 5630544.364, 682378.978 5... | 57 |
| 2577 | 2153.0 | POLYGON ((680458.914 5630872.499, 680464.377 5... | 58 |
| 2578 | 568.0 | POLYGON ((681002.891 5632199.397, 681000.911 5... | 59 |
| 2579 | 320.0 | POLYGON ((683485.762 5631278.394, 683479.331 5... | 61 |
2580 rows × 3 columns
streets_eti_prcd
| enclosure_tess_index | is_first_line | |
|---|---|---|
| 0 | 0 | 1 |
| 1 | 1 | 1 |
| 2 | 2 | 0 |
| 3 | 3 | 1 |
| 4 | 4 | 1 |
| ... | ... | ... |
| 2575 | 2575 | 1 |
| 2576 | 2576 | 0 |
| 2577 | 2577 | 0 |
| 2578 | 2578 | 0 |
| 2579 | 2579 | 0 |
2580 rows × 2 columns
# добавляем к полигонам геодатафрейма enclosed_tess признак is_first_line
result = enclosed_tess.merge(streets_eti_prcd[['is_first_line']], left_index=True, right_index=True)
# визуализируем границы enclosed tesselations и enclosures (красная линия), подсвечены полигоны, которые примыкают к улицам с типом 1
ax = result.plot(column='is_first_line', edgecolor='white', linewidth=0.2, figsize=(10, 10))
buildings.plot(ax=ax, color='white', alpha=.5)
enclosures.plot(ax=ax, edgecolor='red', linewidth=2, alpha=.2)
ax.set_xlim(681500, 682500)
ax.set_ylim(5631000, 5632000)
ax.set_axis_off()
gdf_start_point = result.to_crs('EPSG:4326')
start_point = gdf_start_point.geometry.unary_union.centroid.coords.xy
location = float(start_point[1][0]), start_point[0][0]
m = folium.Map(location=location, zoom_start=14, overlay=True, show=False)
folium.TileLayer('', name='Отключить слои', attr='blank').add_to(m)
df_streets.explore(m=m,
column='highway',
name='Дорожная сеть',
overlay=True,
legend=False,
show=True
)
result.explore(m=m,
column='is_first_line',
name = 'Первая линия домов или нет',
overlay=True,
legend=False,
show=True
)
folium.LayerControl(collapse=False).add_to(m)
<folium.map.LayerControl at 0x2c1b1320e60>
m